function res = imo_solve(P,params)
    % OBJECTIVE : Solve the equilibrium, given starting value P and params.
    
  function Pind = Priceind_imo_2009(x_x,params)
    % Objective: Given the cutoffs, what are the price indices?
    K = params.K;
    sigma = params.sigma;
    phi = params.phi;
    x = params.x;
    Tau = params.Tau;
    t = params.t;
    
    for j=1:K       % Calculate the price index for country j
      for k=1:K     % Find country k's contribution to P in country j
            x_cut_x = x(x>x_x(k,j));        % Productivities above the cutoff
            
            % If no firms are over the cutoff, the integral is 0. We set a value close to zero,
            % so that we get a value for P even if no firms present from
            % any source countries.
            if length(x_cut_x)==0 tmp(k)=1e-5;
            else                
                markup = sigma/(sigma-1);
                price = markup.*(Tau(k,j)./x_cut_x+params.t(k,j)); 
                share = length(x_cut_x)/length(x);                     % Pct. of firms above cutoff
                integral = share * mean(price.^((1-sigma)*(1-phi)));
                tmp(k) = params.potential_entrants(k)*integral;        
            end
      end
      Pind(j,:) = sum(tmp)^(1/(1-sigma));
      
    end
  end

  function Pind = Priceind_closed_form_imo_2009(x_x,params)
    % Objective: Given the cutoffs, what are the price indices?
    % This version is using Pareto and iceberg costs only (fast)
    K = params.K;
    sigma = params.sigma;
    gamma = params.gamma;
    Y = params.potential_entrants;
    mu2 = (sigma/(sigma-1))^(1-sigma)*gamma/(gamma-(sigma-1));
    %x = params.x;
    Tau = params.Tau;
    for n=1:K       % Calculate the price index for country n
      for i=1:K     % Find country i's contribution to P in country n
            tmp(i) = Y(i).*Tau(i,n).^(1-sigma).*x_x(i,n).^(sigma-gamma-1);
      end
      Pind(n,:) = (mu2*sum(tmp)).^(1/(1-sigma));
    end
  end


% Find x_x and pi simulaneously

pi = solve_pi_xx(P,params);
Yy = params.L*(1+pi);
x_x = cutoff_imo_x(P,Yy,params);

P2 = Priceind_imo_2009(x_x,params);
%P2 = Priceind_closed_form_imo_2009(x_x,params);

% fsolve. We have solved the system when res is 0.
res = P2 - P;

end
